arXiv: 1502.03254v2 [q-fin.PR] 22 Nov 2016 


MASS AT ZERO IN THE UNCORRELATED SABR MODEL AND IMPLIED 

VOLATILITY ASYMPTOTICS 

ARCHIL GULISASHVILI, BLANKA HORVATH, AND ANTOINE JACQUIER 


Abstract. We study the mass at the origin in the uncorrelated SABR stochastic volatility 
model, and derive several tractable expressions, in particular when time becomes small or large. 
As an application-in fact the original motivation for this paper-we derive small-strike expansions 
for the implied volatility when the maturity becomes short or large. These formulae, by definition 
arbitrage free, allow us to quantify the impact of the mass at zero on existing implied volatility 
approximations, and in particular how correct/erroneous these approximations become. 


1. Introduction 

The stochastic alpha, beta, rho (SABR) model introduced by Hagan, Kumar, Lesniewski and 
Woodward in [24[ I26j is now a key ingredient-and has become an industry standard-on interest 
rates markets mm [3 Eg. It is defined by the pair of coupled stochastic differential equations 

d X t = Y t xfdW t , X 0 =x 0 > 0, 

(1-1) d Y t = vY t dZ u Y 0 =y o >0, 

d (Z, W) t = pdt, 

where v > 0, p € (—1,1), /? € (0,1), and W and Z are two correlated Brownian motions on 
a filtered probability space (f2, F, (F t )t>o, P). Its popularity arose from a tractable asymptotic 
expansion of the implied volatility (derived in [24]), and from its ability to capture the observed 
volatility smile; calibration therefore being made easier using the aforementioned expansion. In 
today’s low interest rate and high volatility environment, the implied volatility obtained by this 
very expansion can however yield a negative density function for the price process X in m, 
therefore exhibiting arbitrage. 

This problem of negative density in low interest-rate environments has been directly addressed 
by Hagan et. al [25], Balland and Tran ,7], and Andreasen and Huge [2], who proposed modifi¬ 
cations of the original SABR model. There exist several refinements to the asymptotic formula 
itself: in 05] Obloj fine tunes the leading order, and Paulot [37] provides a second-order term. 
In certain parameter regimes the exact density has been derived for the absolutely continuous 
part (on (0,oo)) of the distribution of X: in the uncorrelated case p = 0, formulae were obtained 
in 0 [28] by applying time-change techniques. The correlated case is much harder, and approxi¬ 
mations have been derived using projection methods in 00, and using geometric tools in [23] . 
Barring computational costs, availability of the distribution of the SABR process is equivalent to 
computing any European prices. This however calls for a computation, not only on the continuous 
part of the distribution (on (0, oo)), but also of its singular part at the origin. Absorbing boundary 
conditions at the origin ensure ensure that the forward rate process A is a true martingale, and 
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the singular part can hence accumulate mass, depending on the starting value of the process, the 
parameter configuration and the time horizon. 

The original asymptotic formula typically loses accuracy for long-dated derivatives, when the 
CEV exponent /? is close to zero, or when the volatility of volatility v is large. The parameter f3 
governs the dynamics of the smile, and small values thereof are usually chosen when the asymp¬ 
totic formula fails, namely on markets where the forward rate is close to zero and for long-dated 
options [71124] . Indeed, it comes as no surprise that the orginal formula-which is an asymptotic 
expansion for small values of v 2 T -breaks down for large maturities, but it is well known that the 
reasons for the inconsistencies of the SABR formula are subtler than that. What we highlight here 
is that the mass at zero can be held accountable for the irregularities in this case as well. While 
standard numerical methods proved reliable when the process remains strictly positive, comput¬ 
ing the probability mass of the SABR model at the origin is a more delicate issue. Due to the 
singularity at the origin, usual regularity assumptions ensuring stability of numerical techniques 
(finite differences or Monte Carlo) are violated at this point, and a rigorous error analysis for these 
methods is not (yet) available. In addition, producing reference values becomes computationally 
intensive for short time scales. 

Since much of the popularity of the SABR model is due to the tractability of its asymptotic 
formula, one should aim at preserving it while taking into account the mass at zero. The parameter 
sets p = 0 or /3 = 0 are the most tractable, and in fact (as observed in m) the only ones where 
certain advantageous regularity properties of the SABR process can be expected. We therefore 
concentrate here on the singular part of the distribution for these regimes, that is, we study 
the probability F(Xt = 0) and provide tractable formulae and asymptotic approximations. The 
relevance of these parameter configurations is emphasized by recent results [5], which suggest a 
so-called ‘mixture’ SABR (a combination of the p = 0 and /? = 0 cases) approach to handle 
negative interest rates in an arbitrage-free way. From a modelling perspective, one may question 
the relevance of an absorbing boundary condition at zero in a financial context, where negative 
rates can actually occur. In fact, from a stochastic analysis perspective, when /? = 0 there is no 
need to impose such a boundary condition. Remarkably however, as pointed out in [5], even in 
market conditions where interest rates become negative, the historical evolution of interest rates 
suggests that their dynamics follow processes whose probability distribution exhibit a singularity 
at the origirQ, which makes the computation of the mass at zero rates relevant for these market 
scenarios as well. 

A further application is a direct approximation of the left wing of the implied volatility smile. In 
order to understand the small-strike behaviour of the SABR smile it is essential to determine the 
probability mass at the origin: asymptotic approximations of the implied volatility are available, 
not only for small and large maturities, but also for extreme strikes. Roger Lee’s celebrated 
Moment Formula (32j-subsequently refined by Benaim and Friz [9] and Gulisashvili [22]-relates the 
behaviour of the implied volatility It{K) for small strike K and maturity T to the behaviour of the 
price process X around the origin. De Marco, Hillairet and Jacquier [12], and later Gulisashvili [20], 
showed that when the underlying distribution has an atom at zero, the small-strike behaviour of 
the implied volatility is solely determined by this mass, irrespective of the distribution of the 
process on (0,oo). We shall numerically confirm this in the (uncorrelated) SABR model, using 
approximations of the probability mass, in agreement with [6,. 

In Section [2] we derive explicit formulae for the mass at zero P(At = 0) in the SABR model 
for finite time as well as for large times in the uncorrelated case. Under this assumption, it is 
possible to decompose the distribution into a CEV component and an independent stochastic time 
change. Such time change techniques have been applied to the SABR model in the uncorrelated 
case in m nu Eg to determine the exact distribution of the absolutely continuous part of the 
distribution on (0,oo). Therefore, our formulae complement these by providing the singular part 
of the distribution (see [23110] for more details about time change techniques in stochastic volatility 
models). In Section 12.21 and Section 12.31 we derive asymptotic expansions for the density of time- 
changed Brownian motion—inspired by the works of Borodin and Salminen [lOj , Gerhold [T8] and 


1 That is, rates ‘stick’ to zero for certain periods of time, see [5] for more details. 
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Matsumoto and Yor [33.— which we use to derive the behaviour of the atom at the origin for 
short and large times. Finally, in Section [31 we use these results to determine the left wing (small 
strikes) of the SABR implied volatility. Using the formulae provided in [l2j[20], we highlight the 
fact that some of the widely used expansions exhibit arbitrage in the left wing, and propose a way 
to regularise them in this arbitrageable region. 


2. Mass at zero in the uncorrelated SABR model 

The price process X in dH) is a martingale m Remark 2]. If we consider X on the state space 
[0,oo), the origin, which can be attained, has to be absorbing [231 Chapter III, Lemma 3.6]. For 
two functions / and g , we shall write f(z) ~ g(z) as z tends to zero whenever lim f(z)/g(z ) = 1. 

z—>o 


2.1. The decomposition formula for the mass. In the case where the correlation coefficient p 
is null, the mass at the origin can be computed semi-explicitly. Conditioning on the path of the 
volatility process Y, the resulting process X satisfies the CEV stochastic differential equation 

dA t = Y t X?dW t , 

starting from Xq = xo, where Y is a deterministic time-dependent volatility coefficient, and 
represents, for fixed w £ H, a realisation of the paths of Y. Consider now the simple CEV 
equation dX t = Xfd Wt starting from Xq, and set 


G t := 


V. 


2(l-/3) 


and 


G t := 


X, 


2(1 — 0 ) 


( 1 -/ 3) 2 ( 1 -/ 3 ) 2 ' 

Then G t = Zjt yi ds , where Z is a Bessel process satisfying the SDE [J 


Subsection 1.1] 


dZ t = jdt + 2 V\Z^\dW t , Z 0 


2(1 -« 
x 0_ 

(1-/3) 2 


By Ito’s formula, the process G solves the same SDE, so that Z = G, and therefore X = X j. y> 2ds - 

It follows that X can be obtained from X using the stochastic time change t > f* Y^ds, namely 
X f = Xjt Y 2 ds . Since this time change is independent of X , one can decompose the mass at zero 
of the SABR model into that of the CEV component at zero and the density of the time change: 

rt 


( 2 . 1 ) 


P(X t = 0 ) = J p( 


X r = 0 


') p (i 


Yrds € dr ) dr. 


where the mass at zero in the CEV model is given by (see [T2] or [SU] Section 6.4.1]) 


( 2 . 2 ) 


(X r = O) 


x r = 0 = 1 — F 


X, 


2(l-/3) 

0 


2(1 — 0 ) ’ 2r(^ — l) 2 


with r, the normalised lower incomplete Gamma function: r(v, z) = F(u) 1 u v 1 e 

Remark 2.1. If /? £ [1/2,1) in (11,11) . the origin is naturally absorbing, and the mass at zero is 
given by (12.21) . When (3 € [0,1/2), the solution to (11.11) is not unique, and a boundary condition at 
the origin has to be imposed. Should one consider the origin to be reflecting, the transition density 
would then become norm preserving, and no mass at the origin would be present. However, it is 
easy to see that there is an arbitrage opportunity if the origin is reflecting. Formula (12.21) carries 
over to the case /3 £ [0,1/2) when the origin is assumed to be absorbing, which we shall always 
consider from now on. This is of course in line with [ 1291 Chapter III, Lemma 3.6], mentioned 
above, which states that the origin has to be absorbing for a non-negative supermartingale. 

Since for each s > 0, Y s is lognormally distributed, we can write 


Yrds £ dr = 


J exp (2uZi""/ 2 )) ds £ df) , 


(2.3) 
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where r := -t-, z\ v ^ 2 ^ ■= Z s — ^os; the density of this functional is given by [TUI Formula 1.10.4] 
^0 


(2.4) 


„2i >Z 


(—'/ 2 ) 


ds € dr ] = 


?3/4 


exp 




8 4v 2 i 




1 


4’ 4v 2 r 


dr, 


where the function m is defined as na page 645]: 
(2.5) 

8z 3/2 r(/r + §)e?v 


m y (n,z) = 




J e zcosh(2 u) yU jyj ^—fj, , 2_2sinh(u) 2 ^ sinh(2u)sin ^^ du, 


and where the Kummer function M reads 


( 2 . 6 ) 


M(a, b, x) = 1 + ^ 


fc=i 


i(a + 1)... (a + k — l)x k 
b(b+1)... (b + k — l)k\ ' 


2.2. Small-time asymptotics. We now study the behaviour of the mass at zero P (X t = 0) as 
time becomes small. The main challenge is to provide a short-time asymptotic formula for the 
density of the time change process, for which standard expansion techniques are not applicable. 
The additive functional arising from the density of an integral over the exponential of Brownian 
motion often appears in the pricing of Asian options and is of interest on its own. This density is 
notoriously difficult to evaluate in small time, due to a highly oscillating factor connected to the 
Hartman-Watson distribution [33] [34] and ED Section 4.6]. These numerical issues are discussed 
in [8j, and Gerhold [18] used saddlepoint methods to provide short-time estimates. Because of the 
time change and the complexity of the Kummer function (in the integrand), small-time asymptotics 
of the mass at zero cannot be estimated directly. Instead, we use an inverse Laplace transform 
approach, inspired by |18] , to provide small-time asymptotic estimates for the density of the time 
change. From (12.41) and (12.51) . we introduce the notation y := 2 v 2 t, and we shall alternate between 
the two notations without ambiguity in order to simplify some of the formulations below. 

Remark 2.2. For w := 1/y, the function m has the form m m (-) = c ro / 0 °° e _ “ 2ro / ro (u)du, for 
some c CT and / ro . One might be tempted to use a standard Laplace method to determine the 
behaviour of as w tends to infinity. However, at the saddlepoint u* = 0, attained at the left 
boundary of the integration domain, all the derivatives of the function / ro -appearing as coefficients 
of the expansion-are null, and the method does not apply. 


We now formulate one of the main results of the paper, which characterises the small-time 
behaviour of the mass at zero in the uncorrelated SABR model. For every r, y > 0, let u v denote 
the largest (positive) solution to the equation 


(2.7) 2yi — 1 + 4 uy + 2 log(z/2)\/u — \/it log(u) = 0, 

2 

with z := jpJy- Clearly, u y depends on r, but we shall omit this dependence in the notation. Set 


My-.= 


!°gK) 

3/2 


3/2 


1 - 

8 U 2 y 


( 2 ‘ 8 ) 

16 Uy'* 8 Uy 

The following theorem, based on m and m , provides a short-time estimate for the mass 
at zero. As showed in the proof, the expansion of the integrand is performed using saddlepoint 
analysis and complex contour deformation. Precise error estimates however require substantial 
additional work and new techniques, which we hope to develop in future publication. The numerics 
performed later in the paper strongly confirm our result. 


Theorem 2.3. In the uncorrelated SABR model, the asymptotic equivalence 


3/ 2 e 5 / 4 

,(X ‘ = °>~?7V^ eXP 


v A t 


exp 


jiogK) 


UyV 


g( r ) 


dr, 


holds 


t tends to zero, where g(r ) = P ^X r = 0 J ex P 


Vo 

4 i/ 2 t 


Theorem 12.31 follows from (12.11) and the following assertion. 
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Proposition 2.4. As y (equivalently t) tends to zero, we have (recall that y = 2 v 2 t) 


f 


Y_ 2 ds € dr 


% /2 e 5 / 4 

r 5 / 4 2 7 / 4 v / u7r 


ex P -77: 


Vo 


exp 


logK) 


1 

M ~2 


UyV 


dr 


16 4u 2 r 

The technical part of the proof relies on the following proposition, proved in Appendix IA.1I 
Proposition 2.5. As y tends to zero, the function m y in (|2.5D satisfies 


M„ 


yiexp(i-/r) 

m y yn , z) ~-—-exp 


2t r 


log K) 


M' 


u v y - 


Remark 2.6. The proof of the proposition uses saddlepoint analysis. The saddlepoint u v is the 
solution to (12.711 . but does not admit a closed-form expression; however, as seen in the proof, it is 
possible to expand it as y tends to zero to obtain 


m v (n,z) = 


Vz\log(y)\ 
2 ypK 


exp < - 


l°g(y ) 2 , |log(y)| 


1 


1 - 


1 


1 , 3/2 


+ O (y 3/2 


% 2 y V ^ / L V 2 y J J 

but numerical computations however show that this estimate is not very accurate. 

2.3. Large-time asymptotics. We now concentrate on the large-time behaviour of the mass at 
zero in the uncorrelated SABR model. From pjOj Formula 1.8.4, page 612], the formula 


exp 


(2 


ds € dr ) = 


7-3/2 


exp 


1 


2u 2 7 


dr 


holds, so that the decomposition (12.11) together with (12.31) imply that (recall that r = 


:= lim F(X t = 0) = —%= 

Roo uV2tt . 


1-T 


x. 


2(1-0) 

0 


2(1 - (3) ’ 2r(/3 - l) 2 


r 3/2 exp 


Vo' 

Vo 

2v 2 i 


dr 


(2.9) 


= 1 - 


Vo 


x, 


2(1-0) 

0 


>\[2/k . 


2(1 - p) ’ 2r(/3 - l) 2 


r 3/2 exp 


Vo 

2v 2 7 


dr. 


When p = 0(= p), the SABR model (11.11) reduces to a Brownian motion on the hyperbolic plane 
(up to a deterministic time change), and a simple computation shows that (12.91) simplifies to 

mi 2 / vxo 

= 1-arctan 


r oo fl-o 

p 7T V Vo 

When (3 ^ 0, the integral in (12.91) does not have a closed-form expression. Expanding the expo¬ 
nential factor for small yo, we can however write, for any n £ N, the 7ith-order approximation 


A — 


1 — F 


1 


2(1-0) 


2(1 “ P) ’ 2r(/3 — l) 2 ) vr*/2^^k\ 


Vo 


1 f Vo 


n 2fe+l / 1 


k —0 


k\v\J 27t V 2u 2 / J 0 

22/o(l ~ /?) _ 

( / I—.:; ) VS/VX fc=0 


Eh(- 

fc =0 
2(1-0) 


dr 


1 - r 


2(1 - p) ’ 2r(/3 - l) 2 

{-l) k (Vo{P~ 1) 2 V r ( fc + 1 + 


2u 2 r , 

r -(fe+3/2) dr 


k\ 


„2 2 ( 1 - 0 ) 

V x 0 


(1 + 2k) 


Note in particular that 

( 2 . 10 ) 


P& } = 


2r 1 


( 1+ 2 ^ 20 ) y Q (\-P) 

{■^ 20 ) 


When r tends to infinity, the integrand clearly converges to zero fast enough. Using the properties 
of Gamma functions in J] Chapter 6], the asymptotic behaviour 


1\ r 1 a exp(—1/r) 

a ’rJ ffc) 


i-r 
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holds as r tends to zero, ensuring that the integral is well defined for all n £ N. Theorem 12.71 
below shows how well (and when) the sequence P^ approximates the mass at zero Poo- Using 
the Taylor formula with Lagrange’s form of the remainder, we obtain 


exp 


I/o 

2v 2 7 


= v ° 


k =0 


(-1 ) n+ \-e( y% \ 


k\ \2 v 2 r J (n + 1)! 


n+1 


2n 2 i 




for some 9 £ (0, j/g/(2^ 2 r)). Therefore, 


(2.11) exp 
For any n > 0, set 

(2.12) b n := 


Vo 

2v 2 J 




k—0 


2 \ fe 

I/O 


k\ \ 2v 2 r 


< 


2 \ n + 1 
Vo 


2j/o(l - /?) 


( 271 =^ 7 ) ^* 0 " 


% 2 (/3^1) 2 


(n+1)! \ 2v 2 r 

( 


T n + l + 


2-2/3 


n!(l + 2 n) 


so that from (12.1111 . (12.121) . and the definitions of Poo and P^\ it follows that, for any n > 0, 


(2.13) 


=£(-!)*&* and 


k—0 


< &n+l- 


Theorem 2.7. The following statements hold for the sequence P( n ) in (12.131) : 

(i) if 2/0 (/3 — l) 2 > v 2 x 2 q 1 ~^ , or Uq(/3 — l) 2 = v 2 ^ 1 ^ and | < j3 < 1, then the se¬ 
quence (Poo^) n >o diverges, and hence cannot be an approximation to the mass at zero Poo/ 

(ii) if 2/0 (/3 — l) 2 < n 2 a;g^ 1-/3 \ or y 2 (ft — l) 2 = u 2 Xq ( ' 1 ~ /3 and 0 < ft < | then 


(2.14) 


= pW + 0 in 1+ 2 - 2(3 exp ( —nlog 


2 2(l-/3) 


v x, 


0 


Uo(P - l ) 2 


as n tends to infinity. 


Remark 2.8. Remember that xo denotes the initial value of the stock price or interest rate. For 
all practical and sensible values of the parameters, condition (ii) in the theorem is always in force. 


Proof. From (12.121) . Stirling’s formula for the Gamma function yields, as k tends to infinity, 

k 

llnf'i — R\ - . R / 7 /^ (R — 1 \ 

(2.15) b k 


2/o(l - /3) 


(271=37) u ^ a 


x-t 




,2 2(l-/3) 


From (12.131) . if the conditions of Theorem 12.71' i') hold, the general term of the series X3fc°=o(— 
does not tend to zero, and the sequence Pl2 diverges. If the conditions of Theorem 12. 7f ii) hold, 
then (12.131) and (12.151) imply (12.141) , which completes the proof of Theorem 12.71 □ 


For practical purposes, depending on the conditions for convergence of the sequence (P^)n>o in 
Theorem l2.71 it may or may not be useful to use directly the integral form (12.91) . For (y 0 , v, f}, xq) = 
(0.1,1.0,0.2,0.2) (for which convergence holds), the mass at zero in this case is Poo = 20.833%. 
Using Theorem 12.71 the table below computes the error using the sequence (P»^): 



n = 0 

n = 1 

n = 2 

n = 3 

n = 4 

n = 5 

IP-P^ J I 

6.43E-3 

3.41E-4 

2.13E-05 

1.43E-06 

1.01E-07 

7.29E-09 

Computation time (in seconds) 

6.8E-05 

8.6E-05 

1.3E-4 

1.9E-4 

2.2E-4 

2.6E-4 


and the table below computes the integral (12.91) using the Python scipy toolpack for quadrature; 
the integral is truncated at some arbitrary value R > 0: 
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R = 20 

o 

II 

O 

co 

II 

o 

00 

II 

0$ 

R = 100 

R = 120 

Absolute error 

2.33E-4 

1.07E-4 

6.77E-05 

4.90E-05 

3.81E-05 

3.11E-05 

Computation time (in seconds) 

7.6E-3 

7.9E-3 

8.9E-3 

9.2E-3 

9.6E-3 

9.9E-3 


These results suggest that convergence of the series expansion is extremely fast. In particu¬ 
lar, event the limit (12.101) . with n = 0, yields a very accurate result, which allows for a simple 
interpretation of the impact of each parameter of the model on the large-time mass at the origin. 

Remark 2.9. One could in principle compare these values with Monte Carlo simulations. How¬ 
ever, as far as we are aware, no rate of convergence for such schemes has yet been proved for the 
SABR model, so one may question numbers generated by simulation. In addition, it is known El 
that in the critical region around zero, Monte Carlo methods are prone to a simulation bias. Nev¬ 
ertheless, for a comparison with the above results we included some corresponding values for the 
mass generated by a Monte Carlo algorithm in Section 12.3.11 below. 

2.3.1. Large-time numerics. We provide below some numerics of the large-time mass at zero de¬ 
rived in (12.91) . In particular, we observe the influence of the parameter /3 (Figure[l]) as well as that 
of the starting point xo (Figure [l]) of the uncorrelated SDE (11.11) . As f3 tends to one (from below), 
the mass at zero is diminishing, even for arbitrarily small values of xo . Likewise, as the initial 
value xq increases, the mass at the origin decreases even for /3 = 0. We shall further comment on 
the importance of the mass at the origin in financial modelling in Section [3] below. 


-- j* =0.001 
- - a*, =0.0801 


00 01 0.2 0.3 04 05 06 0.7 08 0.9 



Figure 1. Influence of /3 on the large-time mass at zero in the uncorrelated 
SABR model with (y 0 ,u) = (0.015,0.6) (left) and (ya,v) = (0.1, 1) (right). 



*0 *0 

Figure 2. Influence of the initial value xo on the large-time mass at zero with 
(yo,u) = (0.015,0.6) (left) and (j/o,u) = (0.1,1) (right). This gives a numerical 
interpretation of ‘feeling the boundary’: as we start the diffusion far enough from 
the origin, the mass at zero becomes small. 

With due caution with respect to their validity (Remark 12.91) , we include for comparison a 
sample of values for the mass at zero obtained by Monte Carlo simulations with M = 1000 and 
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M = 2000 paths, and the corresponding computation times, for different time horizons T. The 
results suggest that the ‘large-time’ regime is already achieved for maturities equal to 15 years. An 
explanation for this phenomenon is provided in na Section 4]. As in Section l2T3l above. we used the 
parameters (yo,v, /3 ,xq) = (0.1,1.0,0.2,0.2), for which the exact mass at zero is P ra = 20.833%. 



T = 10 

T = 15 

T = 20 

T = 30 

T = 50 

T = 100 

Monte Carlo mass (M=1000) 

0.1889 

0.2020 

0.2100 

0.2110 

0.2050 

0.2090 

Computation time (in seconds) 

3.222 

3.185 

3.221 

3.179 

3.163 

3.178 



T = 10 

T = 15 

T = 20 

T = 30 

T = 50 

T = 100 

Monte Carlo mass (M=2000) 

0.2100 

0.2075 

0.2050 

0.2100 

0.2065 

0.2175 

Computation time (in seconds) 

6.4437 

6.8386 

6.4805 

7.6308 

6.6587 

6.372 


3. Implied volatility and small-strike expansions 


The implied volatility is the Black-Scholes volatility parameter that allows to match observed 
(or computed) European option prices; it obviously depends on strikes and maturities (see for 
example m for more details). Given a model, the classical route to compute the implied volatil¬ 
ity is (i) to compute the price of the Call (or Put) option, and (ii) to invert the Black-Scholes 
formula. Both steps are numerically demanding, and scarcely provide insights on the behaviour 
of the implied volatility smile. Another route, which has motivated the use of asymptotic meth¬ 
ods in finance, is to obtain closed-form expansions for the smile (for small/large maturities, or 
strikes); however, being asymptotic results, they may lose accuracy when some parameters are not 
small/large enough. The ‘classical’ approximation by Hagan et al [25] is such a formula, which has 
been used extensively by practitioners, despite exhibiting flaws-namely arbitrage-in some regions. 
This anomaly can in principle be fixed if one accounts for the accumulation of mass at zero due 
to the Dirichlet boundary condition. Let us recall a few (model-independent) results regarding 
small-strike asymptotics of the implied volatility. For any strike K > 0 and maturity T > 0, let 
us denote by It(K) the implied volatility. In the presence of strictly positive mass at zero, the 
small-strike tail of the implied volatility satisfies [32] : 


(3.1) 


lim sup 

ki o 


It(K) 



This behaviour was recently refined by De Marco, Hillairet and Jacquier [12], and later by Gulisas- 
hvili [20]. Assuming that P(X T < K) — F(X T = 0) = 0((| log K\~ 3 / 2 ) as K tends to zero, De 
Marco, Hillairet and Jacquier [121 Proposition 3.1] derive the small-strike asymptotic formula 


(3.2) I T (K) = 


2|logA"| J\f 1 (m r ) ( J\f '(m-r)) 2 + 2 ( J\T 1 (m T )) 


T 


VT 


+ 


o 


(i log/x r 3/2 ) 


2^27] log’ A\\ogK\VT 

where m-r := P (Xt = 0) is the mass at the origin, and A f the Gaussian cumulative distribution 
function (an alternative formulation of (13.21) canbe found in 120]]). 


3.1. Comparison with Obloj [36]. Obloj [36] derived an-refined version of the one in .26]- 
implied volatility expansion in the SABR model. However, as we illustrate in Figure 13.11 this 
formula exhibits arbitrage for small strikes. As explained by Gatheral El Proof of Lemma 2.2], 
the density of the log price log(X) (or log forward rate) can be expressed directly in terms of the 
implied volatility, and negative densities obviously yield arbitrage opportunities. In Figure 13.11 
we visually quantify how ‘wrong’ Hagan’s expansion is for small strikes in the presence of a mass 
at the origin. We plot k H > Ir(e k )\/T/\k\, which, from (13.21) has to be bounded by y/2 in order 
to avoid arbitrage, and compare it to the first and second order of (13.21) using (12.91) to compute 
the (large-time) mass at zero. We consider two parameter sets, one for which the large-time mass 
is small, and one for yielding a large mass at the origin. As the mass becomes small, Hagan’s (or 
Obloj’s) approximation becomes more accurate. This holds in particular as the parameter {3 gets 
close to one, as indicated in Section f2. 3. II above. In the limit as (3 = 1, the mass becomes null. 
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Figure 3. Density (right) of the log process log(X) obtained from the implied 
volatility expansion [35] (left) with (v,/3,p,xo,yo,T) = (0,1,0.6,0.05,0.5,1.2). 
The mass at zero, computed using (EH) is equal to 4.5%. 



log strikes log strikes 


Figure 4. The black line marks the level \/2. The parame¬ 
ters are (v, /?, p, xq, Do, T) = (0.3,0,0,0.35,0.05,10) for the left plot, and 
(u, /3, p, xq, 2 / 0 , T) = (0.6,0.6,0,0.08,0.015,10) for the right graph. Obloj’s im¬ 
plied volatility expansion clearly violates this upper bound in both cases. The 
large-time mass is equal to 28.3% for the left plot and 3.1% for the right one. 


3.2. Comparison with Antonov-Konikov-Spector [5]. In the uncorrelated case p = 0, Antonov, 
Konikov and Spector [6] derived the double integral formula for the price of a Call option: 


E(X t -K)+ = (X 0 -K)+ + 


2\JXqK j f s+ sin (r}ip{s)) . 2 


sinh(s) 


G(u 2 T, s)ds + sin(?77r) 


exp(-#(s)) 2 


' s + 


sinh(s) 


G(u 2 T, s)ds 


where rj := l/|2(/3- 1|, q := q 0 ~ s± := arcsinh \^\q±q 0 \j, 


ip(s) := 2 arctan < 


/sinh(s) 2 — sinh(s_) 2 
sinh(s+) 2 — sinh(s) 2 


and ip(s) '■= 2arctanh< 


/sinh(s) 2 — sinh(s+) 2 
sinh(s) 2 — sinh(s_) 2 


The function G is defined as the integral 


G(t,s) 


2 exp (— 1/8) 

t 3 / 2 ^ 


cosh(n) — cosh(s) exp 



dzi. 


In Figure 13.21 we compare the smile obtained from the Antonov-Konikov-Spector’s formula 
(computing the double integral and numerically inverting the Black-Scholes formula) and the 
closed-form tail formula (13.21) using the large-maturity mass at zero computed from (12.91) . Follow¬ 
ing 0, we consider the following set of parameters: (v,/3,p,Xo,yo,T) = (0.8,0.1,0,0.1,0.15,20), 
for which the large-maturity mass at zero ESP is equal to 63%. 
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4. Conclusion 

The SABR model is a pillar of mathematical modelling on fixed income desks, but suffers 
from some issues in low interest rate environments, where the process can hit the origin with 
non-zero probability, creating, in most existing approximations (used by practitioners) arbitrage 
opportunities. In this paper, we endeavour to provide accurate estimates for this mass at zero in 
order (i) to quantify the error made by existing approximations, and (ii) to suggest an alternative 
parameterisation of the implied volatility smile for low strikes, ensuring arbitrage opportunities 
do not occur. 


Appendix A. Proofs of Section [2] 


A.l. Proof of Proposition 12.51 Our proof is inspired by [TS], which is based on an inverse 
Laplace transform approach. From uni Page 645], the Laplace transform of the function m y has 
a closed-form representation, namely, whenever /i > —3/2 and z > 0, 

—i /"r(/x+ \ + /0 ^ 

Z ) £u ( p(i + 2v ^) 

where the function 9)1 is related to the Kummer function M function via the identity 
9 7l n ,m(x) = x m+1/2 exp — nf-,2m+l,ij . 

Therefore, we can write, for some R £ R, 


(A.l) m y {n,z) = 


rR+io ° rou + i + ^i) 


2i7T 


(2 z) 2 +v/ “M yj, + — + y/u, 1 + 2 y/u, 2 z ) dit. 


Ir -ioo r(l + 2- v /u) 

Since we wish to determine the behaviour of m y as y (equivalently, t) tends to zero, we need to 
understand the limit of the integrand as u tends to infinity. The following asymptotic relations 
hold uniformly in v, as v = \Ju tends to infinity: 


M ( n + - + v, 1 + 2v, 2z 


e . 


(A.2) T(1 + v) = V2^e~ v v v+1/2 [l + 0{ z;" 1 )] and 

The first one is standard 131] Section 3.5]. As for the second one, the representation (12.611 yields 


1 


m(^+v + h : 1 + 2v,2z'\ =^ 7 fe - 


k =0 


( 2 z_l 

k\ 


where, for k > 0, 


7 k ■= 


(n + v + |) • ■ ■ (n + v + k — |) 
{l + 2v)(2 + 2v)---{k + 2v) ' 


Clearly | 7 fc| < 2 k and 7 fe ~ 2 k as v tends to infinity, and from [TJ Formula 13.6.3], we have 


M ( - + v + n, 1 + 2v, 2z ) < M ( - + v, 1 + 2v, 2z ) = T(1 + v)e z ( - 


I v(z), 
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where again I„ denotes the modified Bessel function of the first kind uni Page 638], so that, 
using (IA.2I) and [T8] Equation 9], we have, uniformly in v, 


M 


+ v + p, 1 + 2v, 2 z^j < r(l + v)e" ^ I v(z) — e‘ (l + O (y x )). 


Therefore the integrand in (IA.1I) reads, as it tends to infinity, 

*(«, V, z) = e uy T ^^J^\ 2z)i+^M ^ + l+^,l + 2VU, 2 z 

e v2 y+ v + z z v +h v v- v - : k 2~ v 


= exp 


v 2 y + av+ ( p ~ z ~ v ) log(v) + 2 + ^ log(z) 


=: exp (ip y (u) + z + - log(z) . 


where a := 1 + log(z) — log(2) € R, and where the function i/j y is defined by 

(A.3) il> y (u) = uy — ^Vu\og(u) + oty/u + \ lo &( u )- 

For y > 0 small enough, the saddlepoint equation d u ip y (u) = 0, or 

2p — 1 + 4 uy + 2(a — 1 )y/u — \/ulog(ii) = 0, 

(namely (12.71) 1 admits a solution u y > 0. This saddlepoint equation can be rewritten as 

log (%) l-a (p- 1/2) 


(A.4) 


y = 


2lly 


Remark A.l. Note that the saddlepoint equation also reads 

= log(uo) _ P _ dp — 2 
y 2yJ2u 0 V2«o 4uo 

where p := \og{z/y/2) and u o := 2 u y , which is reminiscent of that of JTH] . In fact, the saddlepoint 
equation above does not admit a unique solution; in order for the latter to be continuous (as a 
function of y), one should take the largest solution. 

Following jl8j , we deform the integration contour in (IA.1I) around the saddlepoint u y to obtain 

—z rR+ioo r; ^R+ioo r; ru y +ioo 

(A.5) m y {p,z) = — $(u,y,z)du~—^ e^ {u) d u ~—l e^ M du, 

JR— ioo 2l7T JR—ioo ^17T J Uy —ic )0 

as y tends to zero. Let A denote the real integration variable, so that u = u y + iA. Around the 
saddlepoint (A = 0), we have the uniform Taylor series expansions: 


yju = ^T y - 


iA 


A 2 


2\/% 8 uT 


O 


A 3 

5/2 


iA A 2 


r — (2 + log(it y ))iA log(u !/ )A 

V^logu = yfu^\ogu y + - a 1 


8 u 


3/2 


O 


log it = log u y A — T 
(1 + log(rtj / ))A 3 


u y 2u y 


O 


so that 


(A.6) ij)y{u) = u v y + y/u^ (^a- 


2y/Uy 

log (%)\ log (Uy) p\og(u y ) 


5/2 


- My A 2 + O 


A 3 (l +log(u y )) 


5/2 


where the coefficients in front of A cancel out from the saddlepoint equation (12.71) . and where 


My,= 


log(M y ) a 1 — 2p 


16 u 


3/2 


8 u- 


3/2 


8 U l 


as defined in Proposition 12.51 By bootstrapping (see Section [A. 1.1 1 for details), the expansion 


(A.7) 


log(l/) 2 
4 y 2 


\ _ 2 log log(1 /y) log(^ 2 ) q 


log (y) 


log (y) 


log (y) 










































12 


ARCHIL GULISASHVILI, BLANKA HORVATH, AND ANTOINE JACQUIER 


holds for the saddlepoint as y tends to zero, and implies 


(A.8) 


y 

M y = r~l~Y2 

log (y) 2 


1 + 0 


{ log I log(y) I 
V log (y) ) 


Since (IA.5I) can be rewritten as 

nUy -\-lh 




\/z 

2+K Ju v — ih 


e^(“)dr 


27r 


exp 


UyV 


a — 


log (u y )\ log(%) , /Gog (Uy) 


e- M « A2 dA, 


' —h 


we need to determine an estimate for the last integral on the compact interval [—h, h]. As explained 
below-where the tail integrals are taken into account-the choice h := log (y) 2 /y 3 ^ 2 is in fact the 
right one, and it is then easy to show that, as y tends to zero, 


/: 


e -M y \ 2 dx = V^l l0 S(2/)l 


,,3/2 


o 


(e- logfe)2 i/- 3/2 ) , 


m y (n,z) 


■ exp 


which then implies 
27T 

2n 


= —— exp 


u v y 

1 

2 " 


a — 


log(u y )\ log(n y ) n\og{uy) 


V^|log(j/)| 


1/3/2 


log(u y ) 


V - o - u vV + 


iAI iog(y)l 

?/3/2 


eXp > 2 


(A.9) 


‘lypK 

z^\log(y)\ 

2 v^r 


|log(y)| 1(^-1) i , ^ 

, i3/2 Uy exp {-u v y + y/u^) 


exp 


log(y) 2 


I log(y)| 


2 _# * 


l-jlog 


log(y) 2 

4y 2 


— r + O (y 

7/2 \ 


4 y 2 y 

where we used the saddlepoint equation (IA.4I) in the fourth line. 

It now remains to prove that one can indeed neglect the tails of the integration domain, where 
3 (m) = A > h. The analysis of this is similar to that of jTS] Section 3], and we only outline here 
the main arguments. First, specify a choice h := log {y) 2 /y 3 ^ 2 of integration bounds accounting for 
the main contribution to the integral exp(i/) y (u))du, with ip y defined in (IA.1I) and where 

u y denotes the saddlepoint in (IA.4I) . By symmetry, it is clearly sufficient to consider only one side 
of the tails, and we shall therefore focus on the positive one e^ y ^du. The analysis is then 

split into studying the inner tail h < A < exp(log(l/i) 2 /4) and the outer tail A > exp(log(l/t) 2 /4). 
Similarly to os Equation (10)], the estimate 

j£ + £ e^du ~ 2 exp j u y t + i log)?/) 2 - exp | 

prevails for the outer tail. For any real number B, OS Lemma 1] remains valid for the behaviour 
of the real part of y/u\og(u) + Byfu with respect to |3(u)|, which allows to bound above the 
inner tail by the value of the integrand at A = h of — M y X 2 \\ = h ~ — 4log(y) 2 multiplied by the 
length of the integration path, which is of order e 10 ^ 1 /*) / 4 ; the relative error is therefore of order 
exp(—± log(y) 2 + o(log(y) 2 )). 

The final part of the error analysis in the expansion (IA.9I) follows from analogous estimates 
to pO Table 1], and the total (both tails) error resulting from the completion to Gaussian integral 

2 exp(-4w 2 ) 


_ _ , exp — -oj ,_ ,_ 

\j 2My Jhy/2My V ^ J yj 2My 

The error 0(X 3 /u^ 2 ) from the local expansion (IA.6I) for is of order 
(A.10) O (log (y) 2 ^y ), 


dw 


= exp 


uj—h\/Mi 


~2 log(t) 2 + o(log(t)) 
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which is immediate from bootstrapping (IA. 8 I) . (IA.7I) for M y and u y , and from the choice of h , 


A 3 < c }og{u y ) 1 log {yf 


5/2 


3/2 


h y 


3/2 


C\og(y) 2 ^/y. 


U,y Uiy 

Hence the total relative error is dominated by the error (1A.10I) from the local expansion if M y is 
not expanded, and by the relative error (I A.8f) of M y if one consider its bootstrapping expansion. 

A. 1.1. Justification of the expansion for M y . We define 


M y := 


log (u y ) a 


16 u' 


3/2 


8 u 


3/2 ' 


The term (1 — 2/z)/(8 u 2 ) in the definition (12.81) of M y is of higher order, so that we can work with 
the simpler expression M y in the bootstrapping expansion and the error analysis instead of M y . 
With a = p + 1 — i log(2) and u = u y / 2, the approximation of the saddlepoint simplifies to 


My = 


log (u v ) a log (u y ) p+1 log( 2 ) 


16 uj 2 


8 uj 2 


16 uj 2 


8 uT 


1 ( \/ 21 og (u) y/2(p + 1 - log( 2 ))^ 


16 uj 2 4 I 16m3/2 


8 u 3 / 2 


Thus M y is up to constants of the same form as [THJ Equation (12)]. By bootstrapping, 


M„ 


log (y) 2 

Indeed, the saddlepoint equation El) 

log( 2 %) 


1 + 0 


y = 


log I log(y)| \ 
log(y) ) 


Ap-2 


2y/2(2uy) y/W^v) 4 ( 2 %)’ 

when setting c(u) = > P = l°g ("Tf) ’ k := Ap — 2 and uq = 2 u y becomes 

y/uo = 2 / _ 1 c(rto), where log(y / uo) = log(c(rto) — log(y)). Hence, bootstrapping as in |l 8 ] yields 


u o = 


1 /log ( 1 / 2 /) k>g(c(u 0 )) 


V 2 \ V2 


P 

V2 


= J_ ( / log(l/y) V 


V2 V2 ' V«o> 

/ log ( 1 /y) \ /log(c(u 0 )) 


V2 ) V V2 J 


V2 


/log(c(w 0 )) p 


y/2 4 y/uf J \ y/2 


V2 4y 'u^ / 


(iog(i /y)f 

1 + 

f ^ ) 

flog (c(u 0 )) p k \ 

, 2 

/Hog(c(u 0 )) p k \ 

to 

1x3 

\}og(l/y) J 

\ y/2 V2 4 v^J 

' (log(l/y )) 2 

{ y/2 y/2 Ay/uoJ 


Now expanding around log(l/y), 

log(c(u 0 )) log (log(l/y)) log( 2 ) 


log (c(u 0 ) 




2y/2uo J 


V2 

and using the fact that both 


V2 


2y/2 


V^logfl /y) 


2V2 


log (c{u) 


- p + 


2-\/2uq 1 


and 


f log (c(u 0 )) - P 


log (y) \ 4 yfu \/ 21 og(y) 

are of order o (l/log(y)), we obtain, by collecting terms, 


log (y) 2 V V2 


+ 


Ay/uo 



2 Uy 

Similarly, 

3/2 1 

u 0 = 3 

yi 

- log(y) 

V2 


log (y) 2 
2 y 2 


1 - 


2 l°g(— log(y)) , 2p + log( 2 ) 


log(y) 


log(c(u)) 

V2 


1 3 


+ 


y/2 Ayfu 


log(y) 


^ l°g(y ) 3 


log(y) 


2 Viy 3 


1 - 


log(-log(y)) 2p + log( 2 ) 


log (y) 


2 log( 2 /) 


log(y) J. 












































































14 


ARCHIL GULISASHVILI, BLANKA HORVATH, AND ANTOINE JACQUIER 


hence u^ 2 ~ (log (1/y)) 2 /(8y 2 ); further, 

log(uo) = -2 (log(y) - log(c (w)))-2 log(y) + 2 log(- log(y)) - log(2) 


2 log (e(n)-p +i ^=) 
log(y) 


so that, by bootstrapping we also recover the form of m Equation (13)], at u = \u y \ 


1 

j2log(u) 

v / 2(p+ 1 - log(2)) 

y 3 

\ , ^ flog(-log(y))Y 

4 

16m 3 / 2 

8 u 3 / 2 

2 log (y) 2 

L + l iog(») )\ 
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